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The  Effect  of  Spray  Initial  Conditions  on  Heat  Release  and 
Emissions  in  LDI  CFD  Calculations 
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NASA  Glenn  Research  Center,  Cleveland,  Ohio,  44135 

Farhad  Davoudzadeh 

Air  Force  Research  Laboratory,  Edwards  AFB,  CA  93524-7680 


The  mass,  velocity  distribution,  droplet  size  and  distribution  of  liquid  spray  has  a  primary 
effect  on  the  combustion  heat  release  process.  This  heat  release  process  then  affects 
emissions  like  Nitrogen  Oxides  (NOx)  and  Carbon  Monoxide  (CO).  Computational  Fluid 
Dynamics  gives  the  engineer  insight  into  these  processes,  but  various  setup  options  exist 
(number  of  droplet  groups,  intial  droplet  temperature)  for  spray  initial  conditions.  This 
paper  studies  these  spray  initial  condition  options  using  the  National  Combustion  Code 
(NCC)  on  a  single  swirler  Lean  Direct  Injection  (LDI)  flame  tube.  Using  laminar  finite  rate 
chemistry,  comparisons  are  made  against  experimental  data  for  velocity  measurements, 
temperature,  and  emissions  (NOx,  CO). 


Introduction 

The  use  of  combustion  Computational  Fluid  Dynamics  (CFD)  in  the  development  of  combustion  technology  has 
been  greatly  facilitated  by  the  advancements  made  during  the  last  decade  in  the  areas  of  combustion  modeling, 
numerical  simulation,  and  computing  platform.  Further  development  of  verification,  validation,  and  uncertainty 
quantification  will  profoundly  impact  the  reliability  and  utility  of  these  modeling  and  simulation  tools.  Under  the 
NASA  Fundamental  Aeronautics  Program,  an  assessment  of  existing  computational  tools  for  emissions  and  flow 
field  is  being  carried  out.  As  a  first  step,  the  present  effort  aims  at  establishing  the  baseline  for  prediction  methods 
and  experimental  data  for  Lean  Direct  Injection  (LDI)1'  2  3  combustion  in  confined,  swirling  flows.  Combustion 
codes  based  on  Reynold  Averaged  Navier-Stokes  (RANS),  Very  Large  Eddy  Simulation  (VLES),  and  traditional 
Large  Eddy  Simulation  (LES)  will  be  used;  the  present  paper  reports  the  preliminary  investigation  using  the 
National  Combustion  Code  (NCC). 

The  National  Combustion  Code  (NCC)  is  a  state  of  the  art  CFD  program  specifically  designed  for  combustion 
processes.  A  short  summary  of  the  features  of  the  NCC  pertaining  to  this  paper  are:  the  use  of  unstructured  grids4, 
massively  parallel  computing  -  with  almost  perfectly  linear  scalability  ,  a  dynamic  wall  function  with  the  effect  of 
adverse  pressure  gradient7,  low  Reynolds  number  wall  treatment8,  and  a  cubic  non-linear  k-epsilon  turbulence 
model9  10,  lagrangian  liquid  phase  spray  model11,  and  stiff  laminar  chemistry  integration.  Recently,  viscous  low- 
speed  preconditioning1213  has  been  added  to  improve  the  low-speed  convergence  of  the  NCC  in  viscous  regions,  and 
the  ability  to  handle  multiple  sets  of  periodic  boundary  conditions  has  also  been  added.  The  combination  of  these 
features  is  usually  not  available  in  other  CFD  codes  and  gives  the  NCC  an  advantage  when  computing  recirculating, 
turbulent,  reacting,  spray  flows.  Previously,  the  NCC  has  undergone  extensive  validation  studies  for  simple 
flows14'15,  complex  flows16,  NOx  emissions  prediction  performance17,  and  traditional  gas  turbine  combustor  / 
injectors18. 

This  paper  extends  the  LDI  combustion  CFD  analysis  of  Davoudzadeh.19  2"  Where  previous  work  only  looked  at 
non-reacting  flows,  this  paper  will  compare  the  NCC  using  a  steady-state  RANS  approach  against  experimental  data 
for  a  confined,  reacting  spray  flow.  This  paper  will  show  the  sensitivity  to  the  spray  initial  conditions,  particularly 
the  number  of  droplet  groups  used  to  define  the  distribution,  and  the  initial  temperature. 
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I.  Flow  Modeling 


A.  Geometry  and  Mesh  Generation 

The  single-element  LDI  swirler  configuration  is  illustrated  in  Figure  1.  Each  element  consists  of  an  air 
passage  with  an  upstream  air  swirler  and  a  converging-diverging  venturi  section.  The  fuel  is  injected  through  the 
center  of  swirler  and  the  fuel  tip  is  at  the  throat  of  the  venture.  The  air  swirlers  have  helical,  axial  vanes  with 

O 

downstream  vane  angles  of  60  .  There  are  six  vanes  with  an  inside  diameter  of  9.3  mm  and  an  outside  diameter  of 
22.1  mm.  The  air  then  dumps  into  a  50.8  by  50.8  mm  combustion  section. 

Since  the  NCC  allows  for  unstructured  elements,  the  grid  may  be  composed  of  any  type  and  mix  of  three- 
dimensional  elements.  However,  hexahedral  elements  were  chosen  because  they  are  more  efficient  at  filling  a 
volume  with  a  smaller  number  of  elements  compared  to  an  all  tetrahedral  grid.  Hexahedral  elements  also  allow  a 
better  calculation  of  the  normal  derivatives  that  are  crucial  for  accurate  boundary  layer  resolution.  Approximately 
850,000  elements  were  used.  The  “Gridgen”  mesh  generation  software  was  used  to  create  all  the  grids  used  in  the 
numerical  simulation  reported  in  this  paper. 

B.  Air  Flow  Conditions 

An  inlet  boundary  condition  is  used,  with  the  air  flow  speed  (20.14  m/s)  normal  to  the  inlet  face,  density  (1.19 
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Kg/m  ),  turbulence  intensity  level  (approximately  5%),  turbulence  mixing  length  (.15m),  and  static  temperature 
(294.28  °K).  The  inlet  pressure  in  these  boundary  conditions  are  calculated,  not  specified  (because  they  are 
subsonic).  The  exit  boundary  condition  for  single-element  specifies  static  pressure.  The  operating  pressure  of  the 
combustor  is  approximately  1  atm,  while  the  equivalence  ratio  is  0.75.  The  measured  pressure  drop  (as  a  percentage 
of  P3)  during  the  experiments  was  measured  at  4  %. 

C.  Liquid  Spray  Conditions 

The  spray  evaporation  process  was  modeled  as  a  lagrangian  particle,  where  each  particle  represents  a  group  of 
actual  spray  droplets.  The  droplets  are  represented  as  a  group  rather  than  an  individual  droplet  because  calculating 
each  droplet  is  impractical,  as  it  is  too  computationally  intense.  From  examining  the  spray  experimental  data,  a 
Sauter  Mean  Diameter  (SMD,  D32)  of  32  microns  was  used  for  all  calculations.  The  hollow  spray  cone  had  a  total 
angle  of  90.0°  and  a  cone  thickness  of  14.0°.  A  mass  flow  of  4.15E-04  kg/s  and  an  initial  velocity  of  20.0  m/s  was 
prescribed.  Simulations  were  performed  with  three,  nine,  and  eighteen  droplet  groups.  Figure  3  shows  the 
cumulative  dropsize  distribution  (Given  by  Raju11)  used  as  the  spray  initial  condition.  Three  droplet  groups  is  not  an 
adequate  representation  of  the  spray  distribution;  this  was  used  as  starting  point  for  calculations.  Using  nine  and 
eighteen  droplet  groups,  the  drop  size  distribution  is  represented  accurately.  A  well  represented  (nine  droplet  groups 
and  above)  distribution  at  a  32  micron  SMD  has  a  droplet  size  range  from  one  to  ninety-five  microns.  Ninety-six 
streams  and  thirty-two  rays  were  used  for  the  spatial  distribution  for  the  spray.  The  lagrangian  solution  process  uses 
an  “unsteady”  spray  model.  Droplet  groups  are  only  integrated  for  a  fraction  of  their  lifetime  (but  restarted  at  this 
point  for  the  next  iteration),  rather  than  completely  to  steady-state.  A  time  step  of  5.0E-07  was  used  as  the 
integration  time  step  (DTML  in  NCC  terminology),  and  an  integration  time  of  1.0E-04  (DGML).  New  droplet 
groups  were  introduced  every  integration  time.  This  process  is  thought  to  allow  better  load  balancing  on  parallel 
computing  platforms,  possibly  allowing  better  coupling  with  the  CFD  solver. 

D.  Chemistry  Model 

Ideally,  we  would  prefer  to  use  detailed  chemical  kinetic  models.  There  are  two  problems  with  this  approach: 
1)  Jet-A  is  a  fuel  and  not  a  substance,  and  there  are  no  universally  accepted  surrogate  fuel  models  for  Jet-A;  2)  the 
computational  costs  associated  with  these  models  make  them  impractical  when  fine  computational  grids  are  used. 
Originally,  a  single-step,  global  chemistry  model  was  used,  shown  in  Table  1.  This  model  was  based  on  propane 
kinetics21,  which  are  close  to  Jet-A’s  reaction  rates.  The  single-step  model  allowed  an  easier  start  up  in  the  solution 
process,  by  reducing  the  computational  requirements  during  the  ignition  phase.  Single-step  models  do  not  allow 
emissions  calculations,  only  heat  release.  Because  of  this,  a  reduced  ten-step,  twelve  species  model  based  on 
propane  kinetics22,23  was  used,  as  shown  in  Table  2.  The  mechanism  was  developed  by  a  gradual  reduction  of 
reaction  steps  and  species  using  sensitivity  techniques.  The  reduced  mechanism  also  describes  the  formation  of 
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Carbon  Monoxide  and  Nitrogen  Oxide.  However,  only  one  nitrogen-oxide  species,  NO,  has  been  used  in  the 
reduced  mechanism.  NO  in  the  reduced  mechanism  represents  the  whole  family  of  nitrogen  oxides  including  nitric 
oxide  by  Zcldoviclr 1  reactions,  prompt  NO  reactions  by  Fenimore25,  and  nitrogen  oxide  formation  through  nitrous 
oxide. 


E.  NCC  Computations 

Staging  was  used  in  the  solution  process;  cold-flow  calculations  and  initial  combustion  calculations  were 
performed  using  a  single-step  chemistry  model  with  lagrangian  spray  until  a  steady  state  solution  was  obtained.  The 
final  stage  of  CFD  calculations  was  performed  by  switching  from  the  one-step  chemistry  model  to  the  reduced 
chemistry  model;  this  was  done  by  changing  the  input  chemistry -parameters  of  the  code.  It  is  important  to  note  that 
no  turbulence  -  chemistry  interaction  model  was  used  for  this  case.  Based  on  past  experience,  this  is  an  adequate 
engineering  modeling  assumption.  ‘ 

The  NCC  computations  for  reacting  and  non-reacting  flow  were  run  in  general  until  the  flow  residuals  were 
reduced  three  orders  of  magnitude.  The  mass  flow  rates  at  the  boundary  conditions  were  also  monitored  as  a 
convergence  criterion.  Dissipation  (JST  type)  was  set  at  0.0  for  second  order  dissipation  (s2)  and  .1  for  fourth  order 
dissipation  (s4).  26  The  value  of  k2,  the  constant  that  scales  the  second  order  dissipation  gradient  switch,  was  set  at 
0.50.  Setting  the  second  order  dissipation  to  zero  is  absolutely  necessary  to  accurately  resolving  flow  features,  like 
jets.  A  CFL  number  of  1.0  was  used.  A  cubic,  non-linear  k-epsilon  model  with  a  variable  Cmu  coefficient  was  used. 
This  model  was  selected  because  of  the  swirling  flow.  A  dynamic  wall  function  with  pressure  gradient  effects  was 
used  to  model  near  wall  turbulent  flow  effects. 

Computations  were  performed  on  a  variety  of  computer  platforms,  namely  SGI  Columbia  supercomputer  at 
NASA  Ames  and  clusters  (Mac  OSX  and  Linux)  at  NASA  Glenn.  The  Columbia  supercomputer  was  preferred 
because  it  was  considerably  faster  because  of  its  high  speed,  low-latency  interconnect.  This  interconnect  was 
important  because  the  Lagrangian  spray  model  created  a  load  unbalance  in  compute  nodes.  The  high  speed 
interconnect  seemed  to  mitigate  the  load  imbalance.  It  takes  approximately  one  week  to  complete  a  single  LDI 
combustion  case  using  64  processors.  After  the  baseline  run  was  completed,  new  runs  starting  from  the  baseline 
case  took  about  two  days  to  one  week,  depending  on  the  number  of  spray  particles  used. 


II.  Experimental  Data 

The  experimental  data  is  provided  by  Jun  Cai,  S.-M.  Jeng,  and  R.  Tacina27,  and  by  Yongqiang  Fu  and  San-Mou 
Jeng.~s  Velocity  measurements  were  taken  with  a  two-component  Laser  Doppler  Velociometry  (LDV)  system, 
temperature  measurements  were  taken  with  thermocouples,  and  emissions  data  was  gathered  via  an  isokinetic  probe 
and  gas  analyzer.  Quartz  makes  up  the  combustion  section.  The  combustor  experiments  have  significant  convective 
and  radiative  heat  losses.  The  temperature  measurements  reported  are  not  corrected  to  adiabatic  conditions.  The 
heat  transfer  losses  will  cause  a  200  to  300  °K  temperature  difference  between  the  experimental  data  and  the  NCC 
computations.  Experimental  droplet  measurements  are  collected  with  a  Phase  Doppler  Particle  Analyzer  (PDPA). 


III.  Results  and  Discussion 


A.  Process  Overview 

Table  3  shows  the  CFD  calculations  performed  to  date.  In  Figure  4,  streamlines  visualize  the  flow  as  it  passes 
through  the  helical  swirler  and  into  the  combustion  chamber.  Note  the  strength  of  the  recirculation  zones  by  the 
track  of  the  streamlines.  Contours  of  axial  velocity  are  also  shown.  Blue  indicates  negative  axial  velocity  (and  the 
recirculation  zones)  and  red  indicates  positive  axial  velocity.  The  recirculation  zone  size  and  strength  affects 
combustion  stability  and  emissions.  Figure  5  shows  the  vaporization  process.  The  liquid  droplets  (small  spheres) 
vaporize  before  combustion  occurs.  The  fuel  vapor  is  shown  by  the  green  iso-surface.  This  fuel  vapor  mixes  further 


'  Results  will  challenge  this  assumption.  However,  the  cost  of  running  a  complex  turbulence-chemistry  interaction 
model  is  very  prohibitive,  even  for  supercomputers  like  NASA’s  Columbia. 
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with  air  and  then  burns.  Contours  of  gas  temperature  along  a  centerline  cut  visualize  this  combustion  process.  Blue 
indicates  room  temperature  and  red  is  at  the  adiabatic  flame  temperature.  Note  how  the  flame  shape  is  bound  by  the 
fuel  vapor  region. 


B.  Flow  Comparision 

Figure  6  shows  a  comparison  of  axial  velocity  at  the  centerline.  All  of  the  simulations  show  general  agreement 
with  the  experimental  data,  but  appear  out  of  phase.  The  experimental  data  also  shows  a  dip  in  axial  velocity  from  0 
-10  mm  from  the  combustor  face.  The  calculations  do  not  produce  this  dip  in  axial  velocity.  Figure  7  gives  some 
hints  to  why  calculations  do  not  agree.  Error  bars  show  the  root  mean  square  (RMS)  values  of  axial  velocity  of  the 
experimental  data  and  the  calculations.  (The  RMS  values  for  the  CFD  calculations  are  assuming  isotropic 
turbulence;  RMS  values  are  the  same  for  each  velocity  and  are  calculated  from  the  turbulent  kinetic  energy.)  The 
RMS  values  are  very  large  in  the  0-10  mm  range.  Given  this  high  turbulent  intensity,  the  turbulence  model  (and 
the  RANS  assumption)  will  not  be  able  to  reproduce  this  result.  Figure  8  shows  a  spanwise  comparison  of  axial 
velocity.  In  general,  the  calculations  trend  better  with  the  experimental  data  further  downstream.  The  3mm  location 
shows  a  complete  disagreement  between  the  experimental  data  and  CFD  calculations.  This  disagreement  was  also 
show  by  the  centerline  comparison.  Both  nine  droplet  group  simulations  match  the  data  better  at  the  3mm  location, 
but  the  reason  is  not  understood.  Figure  9  shows  the  comparison  with  radial  velocity  for  spanwise  rakes.  The  data 
is  reasonable  matched  in  all  locations.  The  simulations  performed  with  more  droplet  groups  tend  to  match  the  data 
better  in  the  3  and  9mm  locations.  Figure  10  shows  a  comparison  of  turbulent  kinetic  energy  along  the  centerline. 
(The  turbulent  kinetic  energy  for  the  experimental  data  was  created  by  duplicating  the  RMS  value  for  the  radial 
velocity  to  get  the  three  values  needed  to  calculate  the  turbulent  kinetic  energy.  The  RMS  of  axial  velocity  is  the 
other  component.)  A  very  interesting  result  occurs  -  changing  the  number  of  droplet  groups  drastically  changes  the 
turbulent  kinetic  energy  produced.  The  nine  droplet  group  simulations  appears  to  qualitatively  match  the 
experimental  data.  Why  this  exactly  occurs  is  unknown.  * 

C.  Temperature  Comparison 

As  a  reminder,  experimental  temperature  is  not  corrected  due  to  heat  losses.  Comparing  temperature  along  the 
center.  Figure  11,  as  great  deal  of  variation  is  produced.  The  one-step  model  shows  a  fast  reaction  rate,  while  the 
ten-step  kinetic  model  shows  that  process  is  a  function  of  the  number  of  droplet  groups  used.  While  the  nine  and 
eighteen  droplet  group  distributions  do  not  have  a  difference  in  droplet  distribution,  the  eighteen  droplet  group 
computation  shows  a  much  faster  reaction  rate.  In  general,  the  CFD  calculations  at  least  qualitatively  match  the 
shape  of  the  experimental  centerline  temperature.  The  more  droplet  groups  used,  the  better  the  agreement.  Using  a 
lower  droplet  temperature  (300°K)  also  seems  to  affect  the  solution  positively.  The  spanwise  rakes.  Figure  12  and 
Figure  13,  also  compare  the  CFD  calculations  against  experimental  temperature.  Generally,  the  comparison  agree 
more  downstream.  Near  the  injector,  there  is  a  considerable  difference.  At  the  5mm  rake,  only  the  calculations  with 
large  numbers  of  droplets  agree  qualitatively  with  the  experimental  data.  However,  the  experimental  data  shows  a 
strong  peak  temperature  at  the  center  of  the  rake.  None  of  the  CFD  calculations  reproduced  this  shape  accurately. 
At  the  10mm  rake,  the  agreement  is  better,  but  only  the  simulations  with  adequate  droplet  distributions  (nine  and 
eighteen  droplet  groups),  give  the  qualitative  shape.  The  shape  of  the  peak  temperature  is  still  not  reproduced 
accurately. 

D.  Emissions  Comparison 

Figure  14  compares  the  amount  of  CO  along  the  centerline.  (Please  note:  the  Y  axis  is  a  log  scale.)  All  of  the 
ten-step  simulations  reproduce  the  general  amounts  of  CO.  The  more  droplet  groups  used,  the  better  the  agreement 
with  the  experimental  data.  Using  a  lower  initial  droplet  temperature  seems  to  positively  affect  the  comparison.  We 
would  like  to  point  out  that  we  have  not  done  any  error  analysis  on  this  data.  Measuring  CO  is  difficult,  because  it 
requires  an  isokinetic  condition  to  prevent  CO  from  oxidizing  to  form  C02.  Spanwise  rakes  of  CO  are  shown  in 


§  While  there  is  no  turbulence  -  chemistry  interaction,  there  is  a  coupling  with  the  heat  release  and  the  turbulence 
equations.  The  heat  release  process  essentially  modifies  the  expansion  and  shear  the  flow  is  experiencing.  This 
shear  then  is  an  input  into  the  turbulence  equations.  We  do  not  claim  that  a  turbulence  -  chemistry  interaction  is  not 
needed. 
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Figure  15.  Once  again,  the  more  droplet  groups  used,  the  better  the  result.  Only  the  simulations  using  nine  and 
eighteen  droplet  groups  qualitatively  agreed  with  the  experimental  data,  particularly  the  shapes.  While  it  appears  the 
CFD  calculations  quantitatively  agree  with  the  experiments,  for  the  reasons  given  above,  this  may  not  be  correct. 
Looking  at  the  NOx  amounts  along  the  centerline.  Figure  16,  we  see  that  the  differences  increase  while  proceeding 
downstream.  Using  more  droplet  groups  does  not  improve  the  CFD  results.  Figure  17  shows  the  spanwise  NOx 
amounts.  Close  to  the  injector  face,  the  CFD  calculations  qualitatively  agree  with  the  experimental  data.  Moving 
downstream  from  the  injector  face,  the  comparison  is  worse.  From  this  comparison,  we  believe  the  ten-step  kinetic 
model  is  not  adequate  for  predicting  amounts  of  NOx . 

E.  Dropsize  and  Mesh  Resolution  -  Sources  of  Error 

Looking  at  Figure  18,  we  compare  computed  Sauter  Mean  Diameter,  after  the  lagrangian  particle  is  injected  into 
the  simulation.  The  CFD  and  experimental  data  are  measured  in  essentially  the  same  way  -  by  probes  to  count  each 
droplet  size.  Details  on  exactly  how  this  is  done  will  be  reported  in  a  later  paper.  While  the  overall  shape  is 
reproduced  by  the  CFD  calculations,  the  dip  in  experimental  droplet  size  is  not  captured  by  the  calculations.  Also, 
the  CFD  calculations  differ  by  more  than  a  factor  of  two.  We  believe  the  32  micron  SMD  given  to  the  distribution  is 
too  small.  For  this  case,  a  70  micron  droplet  size  would  be  more  appropriate.  Figure  19  shows  the  recirculation 
zone,  via  velocity  vectors,  near  the  fuel  injector  tip.  The  strength  of  the  recirculation  zone  is  large  enough  to 
possibly  cause  secondary  droplet  breakup.  For  a  better  comparison  in  the  future,  secondary  droplet  breakup  should 
be  used.  Looking  at  the  gaseous  fuel  mass  fraction  at  the  injector  tip.  Figure  20,  there  is  a  considerable  amount  of 
evaporation  that  occurs.  Computational  grid  also  determines  the  solution  accuracy  near  the  injector  tip.  Figure  21 
shows  the  mesh  spacing  near  the  injector  tip.  Given  the  amount  of  change  that  occurs,  the  mesh  should  be  refined  in 
this  area. 


IV.  Conclusions 

These  CFD  calculations  are  a  step  forward  in  anchoring  CFD  codes,  like  the  NCC,  to  a  swirling,  reacting  spray 
flow.  The  current  results  are  acceptable,  but  can  and  will  be  improved.  The  improvements  that  will  be  made  are: 
using  a  better  reduced  Jet-A  kinetic  mechanism,  using  a  larger  Sauter  Mean  Diameter,  adding  a  secondary  droplet 
breakup  model,  improving  the  mesh  density  near  the  injector  tip,  and  using  a  consistent  turbulence  -  chemistry 
interaction  model. 
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Figure  1.  Single  element  LDI  injector  geometry. 
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Figure  2.  Computional  grid  for  the  single  element  LDI  combustor. 
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Table  1  Single  Step  (Global)  Chemistry  Model. 


Reaction 

A 

(mole  -  cm  -  sec-  K) 

n 

E 

(cal/mole) 

1 

4  C12H23  +  71  02  =>  48  0  02  +  46  H20 

GLO/C12H23  0.10/ 

GL0/02  1.65/ 

8.60E+11 

0.00 

3.00E+4 

Table  2.  Reduced  10  Step,  12  Species  Chemistry  Model  in  Chemkin  Format. 


Reaction 

A 

(mole  -  cm  -  sec  -  K) 

n 

E 

(cal/mole) 

1 

4  C12H23  +  47  02  =>  48  CO  +  46  H20 

GLO/C12H23  0.1  / 

GLO/02  1.6/ 

1.46E+13 

0.00 

3.40E+4 

2 

H2  +  02  <=>  H20  +  O 

3.98E+11 

1.00 

4.80E+4 

3 

H2  +  O  <=>  H  +  OH 

3.00E+14 

0.00 

6.00E+3 

4 

H  +  02  <=>  O  +  OH 

4.00E+14 

0.00 

1.80E+4 

5 

CO  +  OH  <=>  C02  +  H 

1.51E+07 

1.28 

-7.58E+2 

6 

H20  +  02  <=>  20  +  H20 

3.17E+12 

2.00 

1.12E+5 

7 

CO  +  H20  <=>  C02  +  H2 

5.50E+04 

1.28 

-1.00E+3 

8 

N2  +  O  <=>  N  +  NO 

1.00E+14 

0.00 

7.50E+4 

9 

N  +  02  <=>  NO  +0 

6.30E+09 

1.10 

6.28E+3 

10 

N  +  OH  <=>  NO  +  H 

3.80E+13 

0.00 

0.00E+0 
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SMD=32,  dD=0.2,  Dmax=3*SMD=96  micrometers 


Figure  3.  Cumulative  drop  size  distribution  (spray  initial  condition). 


Table  3.  Summary  of  CFD  calculations. 


Chemistry 

Model 

Number  of  Droplet 
Groups 

Droplet  Initial 

Temperature  [°K] 

1  Step 

3 

350 

10  Step 

3 

350 

10  Step 

9 

350 

10  Step 

9 

300 

10  Step 

18 

350 
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Figure  4  Flow  visualization  using  streamlines  and  contours  of  axial  velocity  at  an  axial  mid-plane. 


Figure  5  Flow  visualization  of  the  heat  release  process  using  the  largrangian  droplet  groups, 
contours  of  temperature  at  an  axial  mid-plane,  and  iso-surfaces  of  the  fuel  mass  fraction. 
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Figure  6.  Comparison  of  axial  velocity  versus  axial  distance  at  the  centerline. 
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Figure  7.  Comparison  of  axial  velocity  with  RMS  values  as  error  bars  versus  axial  distance  at  the 
centerline. 
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Figure  10.  Comparison  of  turbulent  kinetic  energy  versus  axial  distance  at  the  centerline. 
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Figure  11. 
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Comparison  of  temperature  versus  axial  distance  at  the  centerline. 
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20mm  Axial  Plane  10mm  Axial  Plane  5mm  Axial  Plane 


Figure  12.  Comparison  of  temperature  versus  spanwise  distance  in  the  axial  Y-Z  plane  with  Z  =  0  mm. 
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Figure  13.  Comparison  of  temperature  versus  spanwise  distance  in  the  axial  Y -Z  plane  with  Z 
mm. 
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Figure  14.  Comparison  of  Carbon  Monoxide  amounts  versus  axial  distance  at  the  centerline. 
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Figure  15.  Comparison  of  Carbon  Monoxide  amounts  versus  spanwise  distance  in  the  axial  Y-Z 
plane  with  Z  =  0  mm. 
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Figure  16.  Comparison  of  Nitrogen  Oxides  amount  versus  axial  distance  at  the  centerline. 
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Figure  17.  Comparison  of  Nitrogen  Oxides  amount  versus  spanwise  distance  in  the  axial  Y -Z  plane 
with  Z  =  0  mm. 
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Figure  18.  Sauter  Mean  Diameter  (SMD,  D32)  for  an  axial  distance  of  5mm  from  the  combustor 
face,  spanwise  rake  along  the  Y  =  Omm. 
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Figure  19.  Velocity  vectors  colored  by  velocity  magnitude  in  an 
axial  slice  at  the  Y=0mm  mid-plane. 
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Figure  21.  Approximate  visualization  of  the  computational  mesh 
in  an  axial  slice  in  at  the  Y=0mm  mid-plane. 


COMPUTATIONAL  METHOD: 

The  computational  work  performed  to  produce  the  numerical  results  presented  in  this  paper  uses  the  National 
Combustion  Code  (NCC),  developed  at  NASA  Glenn  Research  Center  (GRC)  for  comprehensive  modeling  and 
simulation  of  aerospace  combustion  systems. 

The  focus  in  the  development  of  the  integrated  system  of  computer  codes  has  been  to  calculate  the  fluid, 
thermal,  and  chemical  characteristics  of  real-world  combustors  to  an  appropriate  level  of  accuracy  and 
turnaround  time  desired  by  designers  and  analysts.  The  two  foremost  important  obstacles  to  turnaround  time 
have  been  grid  generation  and  single  processing. 

Uses  of  unstructured  or  overset  grids  and  parallel  computing  minimize  the  overall  time  needed  to  achieve  a 
numerical  solution.  Thus  the  main  focus  has  been  to  incorporate  a  numerical  scheme  that  allows  use  of  large 
number  (thousands)  of  processors  in  parallel  to  shorten  the  solution  time  and  to  provide  speed-ups  that  does  not 
deteriorate  with  addition  of  more  processors. 

The  main  flow  solver  for  the  code  used  in  this  work  is  based  on  an  explicit  four-stage  Runge-Kutta  scheme, 
which  is  very  suitable  for  parallelization.  Figure  1  shows  an  example  of  the  speedup  that  has  been  achieved 
with  the  code  on  an  SGI  Origin  2000  [1].  This  3-D  test  case  uses  1.3  million  tetrahedral  elements  for  simulation 
of  a  premixed  hydrogen/air  combustor  [2],  using  the  Intrinsic  Low  Dimensional  Manifold  (ILDM)  kinetics 
module  [3,4].  The  parallel  speedup  metric  is  calculated  by  taking  the  ratio  of  the  time  per  iteration  for  the  serial 
case  versus  the  time  per  iteration  for  the  parallel  case.  The  parallel  efficiency  is  the  ratio  of  the  parallel  speedup 
to  the  number  of  processors  used  in  the  calculation.  In  the  calculations  presented  in  this  paper,  the  authors 
increased  the  number  of  processors  from  200  to  400  and  achieved  approximately  a  factor  of  two  speedup  for  a 
2.5  million  elements  domain. 

To  ease-up  the  grid  generation  task,  the  code  is  designed  to  use  unstructured  meshes.  It  uses  triangular  and/or 
quadrilateral  elements  in  the  2-D  cases,  and  tetrahedrons,  wedges,  pyramids,  and  hexahedrons  in  the  3-D  cases. 
A  combination  of  these  grid  types  can  be  used  to  create  hybrid  grids.  For  example,  to  resolve  the  boundary 
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layer  one  may  choose  to  use  hexahedron  elements  in  the  wall  region  and  transition  out  of  the  boundary  layer  to 
tetrahedron  elements  via  pyramid  elements. 

In  brief,  the  flow  solver  solves  unsteady,  3-D,  compressible  Navier-Stokes  equations.  The  discretization  begins 
by  dividing  the  computational  domain  into  a  large  number  of  elements,  which  can  be  of  mixed  types.  A  central- 
difference  finite-volume  scheme  augmented  with  numerical  dissipation  is  used  to  generate  the  discretized 
equations,  which  are  then  advanced  temporally  by  an  explicit  4-stage  Runge-Kutta  scheme.  For  low  Mach 
number  compressible  flow,  a  pre-conditioning  is  applied  to  the  governing  equations,  and  the  solution  is 
advanced  temporally  by  a  so-called  “dual-time-step”  approach,  in  which  the  Runge-Kutta  scheme  is  used  for 
the  “inner”  iteration.  The  turbulence  model  used  in  the  present  work  is  a  cubic  non-linear  k-epsilon  model  [8] 
with  low  Reynolds  number  wall  integration.  A  description  of  this  solver  and  some  benchmark  test  cases  can  be 
found  in  Refs.  [5-6]. 


Fig  22.  Speedup  curve  for  the  1.3M 
element  test  case  (Quealy,  2002  [1]) 


COMPUTATIONAL  DOMAIN  AND  THE  MESH. 

The  numerical  simulation  is  performed  for  the  whole  geometry  including  the  flow  development  sections  for  the 
air,  six  swirling  air  passages  for  each  module  and  the  rectangular  section  combustion  chambers. 

Several  grid  densities  are  generated  to  consider  the  grid  effects.  For  the  single-element  geometry  three  grids 
with  different  grid  densities  namely  coarse,  medium,  and  fine  are  used.  The  grid  densities  for  the  coarse, 
medium,  and  fine  grids  are  240,384  elements,  624,384  elements,  and  861,823  elements,  respectively.  The  grids 
consist  of  hexahedron-only  elements.  In  contrast  to  grid  topologies  where  the  centerline  becomes  the  axis  of 
singularity  around  which  wedge  type  elements  are  generated,  in  this  all -hexahedron  mesh,  there  is  no  axis  of 
singularity  (see  Figure  3). 

Since  the  computational  code  allows  for  unstructured  elements,  the  hexahedron  elements  can  be  constructed  in 
any  arbitrary  arrangements  relative  to  each  other.  This  allows  some  degrees  of  freedom  for  creating  the  grid 
which  in  turn  accelerates  the  grid  generation  process.  The  grid  generation  part  of  the  computational  work  took 
approximately  one  week  to  complete.  Although  creation  of  the  all -hexahedron  grids  requires  more  effort, 
relative  to  tetrahedron  grids,  but  it  reduces  overall  number  of  elements  required  to  achieve  the  same  solution 
accuracy. 

The  Gridgen  software  is  used  to  create  all  the  grids  used  in  the  numerical  simulation  reported  in  this  paper. 
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•  Hex-only  Mesh 

•  No  Axis  Singularity 


•  Hex-wedge  Mesh 

•  Axis  Singularity 


Figure  23.  Centerline  region  grid  topologies. 
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